4  TWI y HAND

Cálculo para el Partido de La Plata

4.1 Resumen

Este documento evalúa dos métricas topográficas simples para el mapeo de peligro de inundación: el Índice de Humedad Topográfica (TWI) y la Altura Sobre Drenaje Más Cercano (HAND). Utilizando el Partido de La Plata como caso de estudio, desarrollamos ambos índices mediante metodologías consolidadas con datos de elevación FABDEM de acceso libre y las bibliotecas pysheds y xDEM. Los resultados se contrastan con modelado hidrológico oficial disponible de la Facultad de Ingeniería como punto de comparación, evidenciando que tanto el TWI como el HAND representan alternativas viables para contextos con limitaciones de datos donde no existe modelado hidrológico avanzado disponible. Ambos enfoques, respaldados científicamente pero con limitaciones reconocidas, constituyen soluciones prácticas para evaluaciones preliminares de riesgo de inundación a escala municipal cuando se carece de información más detallada.

4.2 Introducción

4.2.1 ¿Qué es el TWI?

El TWI es un índice establecido que combina la pendiente del terreno con el área de drenaje aguas arriba para identificar zonas propensas a la acumulación de agua e inundaciones. Calcula valores distribuidos espacialmente donde números más altos indican mayor potencial de encharcamiento y valores más bajos sugieren condiciones más secas. Es una herramienta ampliamente reconocida en hidrología para modelar condiciones de humedad del paisaje (Atlas 2025).

4.2.2 ¿Qué es HAND?

HAND (Altura Sobre Drenaje Más Cercano) es un descriptor cuantitativo del terreno que representa la diferencia de elevación entre cada pixel en la ladera y el punto más cercano en la red de drenaje hacia donde fluye (Johnson et al. 2019). Este método ha sido ampliamente utilizado para la predicción de extensión de inundaciones porque produce resultados comparables a marcos de modelado más complejos como HEC-RAS, pero con menores requerimientos computacionales (Johnson et al. 2019). HAND ha demostrado ser especialmente efectivo como indicador computacionalmente eficiente de susceptibilidad a inundaciones, requiriendo únicamente datos topográficos como entrada (Watson et al. 2024).

4.2.3 ¿Por qué usar estas métricas?

Tanto el TWI como HAND son métodos establecidos que han demostrado su validez científica en múltiples aplicaciones. El TWI ha sido utilizado por agencias como el Illinois State Water Survey para identificar áreas urbanas con riesgo de inundación (Ballerine 2017), mientras que HAND ha sido implementado por el Centro Nacional del Agua de Estados Unidos para mapeo de inundaciones a escala nacional (Johnson et al. 2019). La principal ventaja de ambas métricas es que proporcionan información valiosa sobre riesgo de inundaciones a un costo extremadamente bajo comparado con estudios hidrológicos detallados. Son gratuitos, rápidos de calcular usando datos topográficos que suelen estar disponibles, y fáciles de interpretar, convirtiéndolos en excelentes herramientas de planificación inicial. Estudios recientes han demostrado que estas metodologías pueden replicar mapas de inundación de alta resolución como indicadores de susceptibilidad a inundaciones (Watson et al. 2024; Li et al. 2023).

4.2.4 Limitaciones

Es fundamental entender que tanto el TWI como HAND son medidas derivadas puramente del terreno y no consideran factores como infraestructura urbana, sistemas de drenaje, vegetación, o patrones climáticos locales. Por tanto, son herramientas que solo proporcionan una noción del riesgo relativo de inundación. No deben utilizarse para tomar decisiones a nivel de parcela específica. Para HAND específicamente, las diferencias en las características del terreno y las incertidumbres asociadas con la estimación óptima de parámetros pueden afectar la precisión de los resultados (Thalakkottukara et al. 2024). Ambas métricas se correlacionan principalmente con flujo superficial y no pueden capturar interacciones complejas con aguas subterráneas o la dinámica completa de sistemas hidrológicos urbanos.

4.2.5 Uso apropiado

Las investigaciones han confirmado que existe correlación entre valores altos de TWI y reportes ciudadanos de inundaciones urbanas menores, validando su utilidad en contextos urbanizados (Kelleher and McPhillips 2020). Similarmente, HAND ha demostrado ser adecuado para el mapeo de inundaciones en regiones con escasez de datos, proporcionando capacidad predictiva comparable para mapear áreas de inundación durante eventos de lluvia extrema (Thalakkottukara et al. 2024). Ambas métricas son especialmente valiosas para gobiernos municipales con recursos limitados como primera aproximación para identificar áreas de riesgo relativo de inundaciones, desarrollar planes de emergencia y priorizar estudios más detallados en zonas críticas. Proporcionan un punto de partida sólido y científicamente respaldado para la gestión del riesgo de inundaciones sin requerir inversión significativa en estudios especializados.

4.3 Herramientas

4.3.1 PySHEDS

“Pysheds es una biblioteca de Python de código abierto diseñada por Matt Bartos para ayudar con el procesamiento de modelos digitales de elevación (DEMs), particularmente para análisis hidrológico. Pysheds realiza muchas de las funciones hidrológicas básicas ofrecidas por software comercial como ArcGIS, incluyendo delineación de cuencas y cálculo de acumulación.” Aquí, utilizamos PySheds para calcular la acumulación de flujo, que se incorpora en nuestro cálculo del TWI.

4.3.2 xDEM

xDEM “fue creado por un grupo de investigadores con experiencia en análisis de datos de elevación para detección de cambios aplicado a glaciología. Hoy en día, su desarrollo es liderado conjuntamente por investigadores en análisis de datos de elevación (incluyendo financiamiento de NASA y SNSF) e ingenieros de CNES (Agencia Espacial Francesa).” Utilizamos xDEM para todo nuestro procesamiento y cálculos de modelos digitales de elevación más allá de la acumulación de flujo.

4.3.3 FABDEM 30m

Fathom, líder de la industria en modelado global de inundaciones, creó FABDEM específicamente para el desarrollo de sus modelos. A diferencia de otros modelos que conservan la altura de construcciones y vegetación, este utiliza técnicas de inteligencia artificial para eliminar dichas interferencias y mostrar únicamente la topografía del suelo (Hawker et al. 2022). Su desarrollo involucró datos de referencia de alta precisión provenientes de doce países con características climáticas y urbanas diversas, lo que garantiza su aplicabilidad en distintos contextos geográficos. Los datos están disponibles para descarga en este enlace. Para nuestros propósitos, FABDEM resulta especialmente valioso dado que las investigaciones han demostrado que la calidad del modelo de elevación constituye el factor más influyente en la precisión del modelado de riesgo de inundaciones. Al eliminar las distorsiones causadas por elementos como edificaciones y árboles, este modelo nos permite obtener cálculos de TWI más precisos y representativos de las condiciones reales del terreno, aspecto crucial para la planificación municipal del riesgo de inundaciones.

4.4 Análisis

4.4.1 Importar datos

En esta sección importamos los datos de elevación FABDEM. Hemos descargado previamente los tiles necesarios para cubrir el área del Partido de La Plata y los importamos usando rioxarray para crear un modelo digital de elevación fusionado que servirá como base para nuestros cálculos hidrológicos.

Mostrar código 
import geopandas as gpd
import matplotlib.pyplot as plt

from pathlib import Path
import xarray as xr
import rioxarray
from rioxarray.merge import merge_arrays
import xdem
import tempfile
import numpy as np
from matplotlib import colors
import leafmap.leafmap as leafmap
from pysheds.grid import Grid
from jenkspy import jenks_breaks

CRS_ARGENTINA = "EPSG:5349"
CRS_WGS84 = "EPSG:4326"

RUTA_BASE = Path("/home/nissim/Documents/dev/fulbright/ciut-riesgo")
RUTA_DATOS = RUTA_BASE / "notebooks/data"
RUTA_PARTIDOS = RUTA_DATOS / "pba_partidos.geojson"

CMAP = "BuPu"

partidos = gpd.read_file(RUTA_PARTIDOS)
partidos = partidos.to_crs(CRS_ARGENTINA)
la_plata = partidos[partidos["fna"] == "Partido de La Plata"]

# Quitar la isla de La Plata - mantener solo el polígono más grande
geometria_principal = la_plata.geometry.iloc[0]
if geometria_principal.geom_type == "MultiPolygon":
    poligono_mayor = max(geometria_principal.geoms, key=lambda p: p.area)
    la_plata = la_plata.copy()
    la_plata.loc[la_plata.index[0], "geometry"] = poligono_mayor

bbox_la_plata_4326 = la_plata.to_crs(CRS_WGS84).total_bounds

rutas_tiles = [
    RUTA_DATOS / "fabdem/S40W060-S30W050_FABDEM_V1-2/S35W058_FABDEM_V1-2.tif",
    RUTA_DATOS / "fabdem/S40W060-S30W050_FABDEM_V1-2/S36W058_FABDEM_V1-2.tif",
    RUTA_DATOS / "fabdem/S40W060-S30W050_FABDEM_V1-2/S35W059_FABDEM_V1-2.tif",
    RUTA_DATOS / "fabdem/S40W060-S30W050_FABDEM_V1-2/S36W059_FABDEM_V1-2.tif",
]

tiles = [rioxarray.open_rasterio(path, chunks=True) for path in rutas_tiles]
dem_fusionado = merge_arrays(tiles)

dem_recortado = dem_fusionado.rio.clip_box(
    minx=bbox_la_plata_4326[0],
    miny=bbox_la_plata_4326[1],
    maxx=bbox_la_plata_4326[2],
    maxy=bbox_la_plata_4326[3],
)

dem_recortado.plot(cmap=CMAP)
ax = plt.gca()
la_plata_wgs84 = la_plata.to_crs(CRS_WGS84)
la_plata_wgs84.plot(
    ax=ax,
    facecolor="none",
    edgecolor="black",
    linewidth=0.5,
    linestyle="--",
    zorder=5,
)

4.4.2 Calcular acumulación de flujo

Calculamos la acumulación de flujo, que es esencial para el cálculo del índice de humedad topográfica, siguiendo el tutorial de pysheds. No calculamos esto usando xarray porque pysheds no es compatible con xarray, pero luego convertiremos estos datos a xarray para nuestro cálculo del TWI. El proceso incluye acondicionar el DEM eliminando pozos y depresiones, calcular direcciones de flujo, y finalmente determinar la acumulación de flujo.

Mostrar código 
with tempfile.NamedTemporaryFile(suffix=".tif", delete=False) as tmp_file:
    ruta_temporal = tmp_file.name


dem_recortado.rio.to_raster(ruta_temporal)
grilla = Grid.from_raster(ruta_temporal)

dem = grilla.read_raster(ruta_temporal)

valor_nodata = dem_recortado.attrs.get("_FillValue", -9999.0)

# Acondicionar DEM
dem_pozos_rellenos = grilla.fill_pits(dem)
dem_inundado = grilla.fill_depressions(dem_pozos_rellenos)
dem_inflado = grilla.resolve_flats(dem_inundado)


dem_inflado_xarray = xr.DataArray(
    dem_inflado,
    coords={"y": dem_recortado.y, "x": dem_recortado.x},
    dims=["y", "x"],
    attrs=dem_recortado.attrs,
).rio.write_crs("EPSG:4326")


mapa_direcciones = (64, 128, 1, 2, 4, 8, 16, 32)


direcciones_flujo = grilla.flowdir(
    dem_inflado, dirmap=mapa_direcciones, nodata_out=np.int32(0)
)

direcciones_flujo_xarray = xr.DataArray(
    direcciones_flujo,
    coords={"y": dem_recortado.y, "x": dem_recortado.x},
    dims=["y", "x"],
    attrs=dem_recortado.attrs,
).rio.write_crs("EPSG:4326")


acumulacion = grilla.accumulation(
    direcciones_flujo, dirmap=mapa_direcciones, nodata_out=np.int32(0)
)


acumulacion_xarray = xr.DataArray(
    acumulacion,
    coords={"y": dem_recortado.y, "x": dem_recortado.x},
    dims=["y", "x"],
    attrs=dem_recortado.attrs,
).rio.write_crs("EPSG:4326")

fig, ax = plt.subplots(figsize=(8, 6))
fig.patch.set_alpha(0)
plt.grid("on", zorder=0)
im = ax.imshow(
    acumulacion,
    extent=grilla.extent,
    zorder=2,
    cmap=CMAP,
    norm=colors.LogNorm(1, acumulacion.max()),
    interpolation="bilinear",
)
plt.colorbar(im, ax=ax, label="Celdas Aguas Arriba")

la_plata_wgs84.plot(
    ax=ax,
    facecolor="none",
    edgecolor="black",
    linewidth=0.5,
    linestyle="--",
    zorder=5,
)
plt.title("Acumulación de Flujo", size=14)
plt.xlabel("Longitud")
plt.ylabel("Latitud")
plt.tight_layout()

4.4.3 Calcular pendiente

Calculamos la pendiente, otra variable necesaria para el cálculo del índice de humedad topográfica, siguiendo el tutorial de xDEM. La pendiente es un componente fundamental de la fórmula del TWI y debe calcularse con precisión en un sistema de coordenadas métricas para obtener resultados confiables.

Mostrar código 
with tempfile.NamedTemporaryFile(suffix=".tif", delete=False) as tmp_file:
    ruta_temporal = tmp_file.name
    
    dem_reproyectado = dem_recortado.rio.reproject(
        CRS_ARGENTINA,
        resolution=30,
    )
    dem_reproyectado.rio.to_raster(ruta_temporal)
    dem = xdem.DEM(ruta_temporal)
    
    atributos = xdem.terrain.get_terrain_attribute(
        dem.data,
        resolution=dem.res,
        attribute=[
            "hillshade",
            "slope", 
            "aspect",
            "curvature",
            "terrain_ruggedness_index",
            "rugosity",
        ],
    )
    
    datos_pendiente = atributos[1]
    
    # Fix coordinate generation - y should go from top to bottom
    coordenadas_y = np.arange(dem.bounds.top, dem.bounds.bottom, -dem.res[1])
    coordenadas_x = np.arange(dem.bounds.left, dem.bounds.right, dem.res[0])
    
    pendiente_xarray = xr.DataArray(
        datos_pendiente,
        coords={"y": coordenadas_y, "x": coordenadas_x},
        dims=["y", "x"],
        attrs={"crs": dem.crs, "units": "degrees", "long_name": "slope"},
    )

    pendiente_xarray.plot(cmap=CMAP)
    ax = plt.gca()
    la_plata.plot(
    ax=ax,
    facecolor="none",
    edgecolor="black",
    linewidth=0.5,
    linestyle="--",
    zorder=5,
    )

4.4.4 Calcular TWI

Combinamos la acumulación de flujo y la pendiente para calcular el índice de humedad topográfica usando la fórmula estándar TWI = ln(α / tan(β)), donde α es la acumulación de flujo y β es la pendiente. Ajustamos los valores extremos que surgen de dividir por pendiente cero para evitar valores infinitos en áreas completamente planas.

Mostrar código 
acumulacion_xarray_reproyectada = acumulacion_xarray.rio.reproject(CRS_ARGENTINA)

# Remuestrear pendiente para coincidir con acumulación
pendiente_remuestreada = pendiente_xarray.rio.reproject(
    acumulacion_xarray_reproyectada.rio.crs,
    shape=acumulacion_xarray_reproyectada.shape,
    transform=acumulacion_xarray_reproyectada.rio.transform(),
)

pendiente_rad = np.radians(pendiente_remuestreada)
datos_twi = np.log(acumulacion_xarray_reproyectada / np.tan(pendiente_rad))

# Reemplazar valores infinitos y valores muy altos
datos_twi = np.where(np.isinf(datos_twi),20, datos_twi)
datos_twi = np.where(datos_twi > 20, 20, datos_twi)

twi_xarray = xr.DataArray(
    datos_twi,
    coords=acumulacion_xarray_reproyectada.coords,
    dims=acumulacion_xarray_reproyectada.dims,
    attrs={
        "crs": acumulacion_xarray_reproyectada.rio.crs,
        "units": "dimensionless",
        "long_name": "Topographic Wetness Index",
        "description": "ln(flow_accumulation / tan(slope + 0.0001))",
    },
)

plt.figure()
twi_xarray.plot(cmap=CMAP)
ax = plt.gca()
la_plata.plot(
    ax=ax,
    facecolor="none",
    edgecolor="black",
    linewidth=0.5,
    linestyle="--",
    zorder=5,
)

Mostrar código 
hand = grilla.compute_hand(
    direcciones_flujo, dem_inflado, acumulacion > 200, nodata_value=np.int32(0)
)

hand_xarray = xr.DataArray(
    hand,
    coords={"y": dem_recortado.y, "x": dem_recortado.x},
    dims=["y", "x"],
    attrs={
        "crs": "EPSG:4326",
        "units": "meters",
        "long_name": "Height Above Nearest Drainage",
        "description": "HAND - Altura Sobre Drenaje Más Cercano"
    },
).rio.write_crs("EPSG:4326")

hand_xarray.plot(cmap=CMAP, figsize=(8, 6))
ax = plt.gca()
la_plata_wgs84.plot(
    ax=ax,
    facecolor="none",
    edgecolor="black",
    linewidth=0.5,
    linestyle="--",
    zorder=5,
)
plt.title("Altura Sobre Drenaje Más Cercano (HAND)", size=14)
plt.tight_layout()

4.4.5 Análisis de distribuciones

Para comprender mejor las características de nuestras métricas, analizamos las distribuciones de valores de TWI y HAND mediante histogramas. Esto nos permite identificar los rangos típicos de valores y la forma de la distribución para cada métrica.

Mostrar código 
# Preparar HAND con transformación de doble raíz cuadrada
hand_xarray_reproyectada = hand_xarray.rio.reproject(CRS_ARGENTINA)
hand_double_sqrt_full = np.sqrt(np.sqrt(hand_xarray_reproyectada))
hand_double_sqrt_mask = hand_double_sqrt_full.rio.clip(la_plata.geometry, la_plata.crs, drop=False, invert=False)
hand_double_sqrt_mask = xr.where(~hand_double_sqrt_mask.isnull(), 1, 0)
hand_double_sqrt_masked = xr.where(hand_double_sqrt_mask == 1, hand_double_sqrt_full, np.nan)
hand_double_sqrt_masked = hand_double_sqrt_masked.rio.write_crs(hand_double_sqrt_full.rio.crs)
hand_double_sqrt_masked = hand_double_sqrt_masked.rio.write_transform(hand_double_sqrt_full.rio.transform())

# Preparar TWI
twi_mask = twi_xarray.rio.clip(la_plata.geometry, la_plata.crs, drop=False, invert=False)
twi_mask = xr.where(~twi_mask.isnull(), 1, 0)
twi_masked = xr.where(twi_mask == 1, twi_xarray, np.nan)
twi_masked = twi_masked.rio.write_crs(twi_xarray.rio.crs)
twi_masked = twi_masked.rio.write_transform(twi_xarray.rio.transform())

# Preparar datos válidos para histogramas
twi_values = twi_masked.values[~np.isnan(twi_masked.values)]
hand_original = hand_xarray_reproyectada.values[~np.isnan(hand_xarray_reproyectada.values)]
hand_double_sqrt = np.sqrt(np.sqrt(hand_original))

# Calcular cortes de cuantiles
twi_quantiles = np.percentile(twi_values, [25, 50, 75])
hand_double_sqrt_quantiles = np.percentile(hand_double_sqrt, [25, 50, 75])

# Obtener colormap para usar los mismos colores que en los mapas
import matplotlib.cm as cm
cmap = cm.get_cmap(CMAP)
cmap_colors = [cmap(i/3) for i in range(4)]  # 4 colores del colormap

# Función para asignar colores por cuantiles usando CMAP
def assign_quantile_colors_cmap(values, quantiles):
    colors = []
    for val in values:
        if val <= quantiles[0]:
            colors.append(cmap_colors[0])  # Q1
        elif val <= quantiles[1]:
            colors.append(cmap_colors[1])  # Q2
        elif val <= quantiles[2]:
            colors.append(cmap_colors[2])  # Q3
        else:
            colors.append(cmap_colors[3])  # Q4
    return colors

# Crear figura con subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

# Histograma TWI con colores por cuantiles
n_bins = 50
counts1, bins1, patches1 = ax1.hist(twi_values, bins=n_bins, alpha=0.8, edgecolor='black')
bin_centers1 = (bins1[:-1] + bins1[1:]) / 2
colors1 = assign_quantile_colors_cmap(bin_centers1, twi_quantiles)
for patch, color in zip(patches1, colors1):
    patch.set_facecolor(color)

ax1.set_xlabel('Valores TWI')
ax1.set_ylabel('Frecuencia')
ax1.set_title('Distribución TWI con Cortes de Cuantiles')
ax1.grid(True, alpha=0.3)

# Agregar líneas de cuantiles
for i, q in enumerate(twi_quantiles):
    ax1.axvline(q, color='black', linestyle='--', linewidth=2, alpha=0.7)
    ax1.text(q, ax1.get_ylim()[1]*0.9, f'Q{i+1}', rotation=90, ha='right')

# Histograma HAND (double sqrt) con colores por cuantiles
counts2, bins2, patches2 = ax2.hist(hand_double_sqrt, bins=n_bins, alpha=0.8, edgecolor='black')
bin_centers2 = (bins2[:-1] + bins2[1:]) / 2
colors2 = assign_quantile_colors_cmap(bin_centers2, hand_double_sqrt_quantiles)
for patch, color in zip(patches2, colors2):
    patch.set_facecolor(color)

ax2.set_xlabel('Valores HAND (√√metros)')
ax2.set_ylabel('Frecuencia')
ax2.set_title('Distribución HAND (Doble Raíz Cuadrada) con Cortes de Cuantiles')
ax2.grid(True, alpha=0.3)

# Agregar líneas de cuantiles
for i, q in enumerate(hand_double_sqrt_quantiles):
    ax2.axvline(q, color='black', linestyle='--', linewidth=2, alpha=0.7)
    ax2.text(q, ax2.get_ylim()[1]*0.9, f'Q{i+1}', rotation=90, ha='right')

# Crear leyenda común usando los colores del CMAP
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor=cmap_colors[0], label='Q1 (0-25%)'),
    Patch(facecolor=cmap_colors[1], label='Q2 (25-50%)'),
    Patch(facecolor=cmap_colors[2], label='Q3 (50-75%)'),
    Patch(facecolor=cmap_colors[3], label='Q4 (75-100%)')
]
fig.legend(handles=legend_elements, loc='upper center', bbox_to_anchor=(0.5, 0.95), ncol=4)

plt.tight_layout()
plt.subplots_adjust(top=0.85)
plt.show()

4.5 Comparativa de modelación FLO-2D con TWI y HAND

Para validar la efectividad de estas métricas topográficas simples, utilizamos datos oficiales de peligro de inundación desarrollados por la Facultad de Ingeniería de la Universidad Nacional de La Plata como parte del Plan de Reducción del Riesgo por Inundaciones en la Región de La Plata (Romanazzi et al. 2019). Estos datos de referencia fueron generados mediante la aplicación del modelo hidrológico-hidráulico bidimensional FLO-2D, que simuló la dinámica de inundación de todas las cuencas del partido de La Plata para distintos escenarios de eventos pluviométricos extremos. La disponibilidad de este modelado detallado proporciona una oportunidad única para evaluar qué tan bien corresponden estos índices topográficos simples con el modelado hidráulico más completo como prueba de concepto.

Para la comparación, presentamos el TWI en su forma continua, permitiendo visualizar la variación completa de valores de humedad topográfica. Para HAND, aplicamos cortes de cuantiles con cuatro clases que dividen el área de estudio en zonas de igual extensión (25% cada una), facilitando la identificación de gradientes de riesgo relativo de inundación. Esta aproximación permite una comparación visual clara entre los diferentes enfoques metodológicos.

Un aspecto importante a considerar es que existen diferencias metodológicas fundamentales entre estos enfoques. El modelado de ingeniería toma en cuenta la infraestructura urbana, especialmente el sistema de calles de la ciudad y elementos del drenaje urbano, mientras que tanto el TWI como HAND se basan únicamente en la topografía del terreno.

Mostrar código 
# Preparar datos para el mapa
ruta_peligro = RUTA_DATOS / "peligro_raster_10m.tif"
peligro_xarray = rioxarray.open_rasterio(ruta_peligro)

# Crear clases de quantiles para TWI usando los breaks calculados
twi_quantile_classes = xr.where(twi_masked <= twi_quantiles[0], 1,
                               xr.where(twi_masked <= twi_quantiles[1], 2,
                                       xr.where(twi_masked <= twi_quantiles[2], 3, 4)))
twi_quantile_classes = xr.where(~twi_masked.isnull(), twi_quantile_classes, np.nan)
twi_quantile_classes = twi_quantile_classes.rio.write_crs(twi_masked.rio.crs)
twi_quantile_classes = twi_quantile_classes.rio.write_transform(twi_masked.rio.transform())

# Crear clases de quantiles para HAND usando los breaks calculados
hand_quantile_classes = xr.where(hand_double_sqrt_masked <= hand_double_sqrt_quantiles[0], 1,
                                xr.where(hand_double_sqrt_masked <= hand_double_sqrt_quantiles[1], 2,
                                        xr.where(hand_double_sqrt_masked <= hand_double_sqrt_quantiles[2], 3, 4)))
hand_quantile_classes = xr.where(~hand_double_sqrt_masked.isnull(), hand_quantile_classes, np.nan)
hand_quantile_classes = hand_quantile_classes.rio.write_crs(hand_double_sqrt_masked.rio.crs)
hand_quantile_classes = hand_quantile_classes.rio.write_transform(hand_double_sqrt_masked.rio.transform())

# Preparar datos de peligrosidad
peligro_2d = peligro_xarray.sel(band=1).astype('float32')
peligro_clipped = peligro_2d.rio.clip(la_plata.geometry, la_plata.crs)

# Definir función para crear mapas consistentes
def create_consistent_map(title, boundary_gdf, bounds=None):
    """Create a map with consistent styling and basemap."""
    fig, ax = plt.subplots(figsize=(8, 6))
    
    if bounds is not None:
        ax.set_xlim(bounds[0], bounds[2])
        ax.set_ylim(bounds[1], bounds[3])
    
    if boundary_gdf is not None:
        boundary_gdf.plot(
            ax=ax,
            facecolor="none",
            edgecolor="black",
            linewidth=0.5,
            linestyle="--",
            zorder=5,
        )
    
    ax.set_title(title, fontsize=14, fontweight="bold", pad=20)
    ax.set_axis_off()
    
    return fig, ax

# Obtener bounds comunes para todos los mapas
common_bounds = la_plata.to_crs(CRS_ARGENTINA).total_bounds

# Mapa 1: TWI
fig1, ax1 = create_consistent_map("Índice de Humedad Topográfica (TWI)", boundary_gdf=la_plata, bounds=common_bounds)

twi_masked_3857 = twi_masked.rio.reproject("EPSG:3857")
twi_masked_3857.plot(ax=ax1, cmap=CMAP, add_colorbar=True, cbar_kwargs={"label": "Índice de Humedad Topográfica"})

plt.tight_layout()
plt.show()

# Mapa 2: HAND
fig2, ax2 = create_consistent_map("Altura Sobre Drenaje Más Cercano (HAND)", boundary_gdf=la_plata, bounds=common_bounds)

hand_quantile_classes_3857 = hand_quantile_classes.rio.reproject("EPSG:3857")
hand_quantile_classes_3857.plot(
    ax=ax2,
    cmap=CMAP + "_r",
    add_colorbar=True,
    cbar_kwargs={
        "label": "HAND (Cuantiles)",
        "ticks": [1.5, 2.5, 3.5, 4.5],
        "format": lambda x, pos: f"Q{int(x)}"
    }
)

plt.tight_layout()
plt.show()

# Mapa 3: Peligrosidad FLO-2D
fig3, ax3 = create_consistent_map("Modelado FLO-2D (Peligrosidad)", boundary_gdf=la_plata, bounds=common_bounds)

peligro_clipped_3857 = peligro_clipped.rio.reproject("EPSG:3857")
peligro_clipped_3857.plot(
    ax=ax3,
    cmap=CMAP,
    add_colorbar=True,
    cbar_kwargs={"label": "Nivel de Peligrosidad"}
)

plt.tight_layout()
plt.show()
(a) Índice de Humedad Topográfica (TWI)
(b) Altura Sobre Drenaje Más Cercano (HAND)
(c) Modelado FLO-2D (Peligrosidad)
Figure 4.1: Comparativa de TWI, HAND y modelado FLO-2D

4.6 Conclusión

Nuestro análisis demuestra que tanto el TWI como HAND corresponden bien al modelado hidrológico desarrollado por la Facultad de Ingeniería Hídrica de la Universidad Nacional de La Plata. A pesar de las diferencias metodológicas entre estos enfoques, la correspondencia espacial general es sólida, proporcionando confianza en el uso de ambas métricas como herramientas de evaluación preliminar para el riesgo de inundaciones. Tanto el TWI como HAND constituyen herramientas valiosas y científicamente respaldadas para evaluaciones preliminares de riesgo de inundación a escala municipal, especialmente en contextos con limitaciones de datos donde no existe modelado hidrológico avanzado disponible. Los resultados son consistentes con estudios internacionales que han validado estas metodologías en ciudades como Katmandú, Nepal (Watson et al. 2024; Li et al. 2023) y en regiones rurales de Estados Unidos (Thalakkottukara et al. 2024), lo que nos proporciona confianza para usar estos datos en evaluaciones iniciales de peligro de inundación.